#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Feb 13 12:31:55 2021

@author: liqingsimac
"""

#6.2.
'''
import numpy as np
xx,yy=np.mgrid[-2:2:41j,-2:2:41j]
zz = -xx/2 -yy/3

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

plt.ion()
fig=plt.figure()
ax=Axes3D(fig)
ax.plot_surface(xx,yy,zz,rstride=5,cstride=4,
                color='c',alpha=0.9)
ax.set_xlim3d(-2.0,2.0)
ax.set_ylim3d(-2.0,2.0)
ax.set_zlim3d(-2.0,2.0)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

#fig.savefig('pic/fig-6-8-1.png')
'''

#6.5.
'''
import numpy as np
import matplotlib.pyplot as plt
x=np.linspace(-2,2,61)
y1=np.exp(x)/2+np.exp(-x)/2
y2=1/y1
fig=plt.figure()
ax1=fig.add_subplot(111)
ax1.plot(x,y1,label='hyperbolic cosine')
ax1.plot(x,y2,label='hyperbolic secant')
ax1.axhline(y=0)
ax1.axvline(x=0)
ax1.legend(loc='best')
'''

## 6.5.1. 交互式显示
'''
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider, Button

def solwave(x,t,c):
    temp=np.cosh(np.sqrt(c)*(x-c*t)/2)
    return c/(2*temp**2)

fig,ax=plt.subplots()
plt.subplots_adjust(left=0.15,bottom=0.30)
plt.xlabel('x')
plt.ylabel('y')
x=np.linspace(-5.0,20.0,1001)
t0=5.0
c0=1.0
line,=plt.plot(x,solwave(x,t0,c0),lw=2,color='blue')
plt.axis([-5,20,0,2])

axcolor='lightgoldenrodyellow'
axtime=plt.axes([0.20,0.15,0.65,0.03],facecolor=axcolor)
axvely=plt.axes([0.20,0.10,0.65,0.03],facecolor=axcolor)

stime=Slider(axtime,'Time',0.0,20.0,valinit=t0)
svely=Slider(axvely,'Vely',0.1,3.0,valinit=c0)

def update(val):
    time=stime.val
    vely=svely.val
    line.set_ydata=solwave(x,time,vely)
    fig.canvas.draw_idle()
    
svely.on_changed(update)
stime.on_changed(update)

resetax=plt.axes([0.75,0.025,0.1,0.04])
button=Button(resetax,'Reset',color=axcolor,
              hovercolor='0.975')
    
def reset(event):
    svely.reset()
    stime.reset()
    
button.on_clicked(reset)

plt.show()
'''

## 6.5.2. 动画
'''
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def solwave(x,t,c=1):
    temp=np.cosh(np.sqrt(c)*(x-c*t)/2)
    return c/(2*temp**2)

fig=plt.figure()
ax=plt.axes(xlim=(-5,20),ylim=(0,0.6))
line, = ax.plot([],[],lw=2)

t=np.linspace(-10,25,91)
x=np.linspace(-5,20.0,101)

def init():
    line.set_data([],[])
    return line,

def animate(i):
    y=solwave(x,t[i])
    line.set_data(x,y)
    return line,

ani=animation.FuncAnimation(fig, animate, init_func=init, 
                    frames=90, interval=30, blit=True)
plt.show()

#To save the animation, use e.g.

ani.save('movie.gif')

#or

writer = animation.FFMpegWriter(
    fps=15, metadata=dict(artist='Me'), bitrate=1800)
ani.save("movie.mp4", writer=writer)

'''

## 6.5.3. 电影任务





##6.7.1. Lissajous曲线，使用 mplot3d
'''
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

theta=np.linspace(0,2*np.pi,101)
a=0.3; m=11; n=9
x=(1+a*np.cos(n*theta))*np.cos(m*theta)
y=(1+a*np.cos(n*theta))*np.sin(m*theta)
z=a*np.sin(n*theta)
 
#plt.ion()
fig=plt.figure()
ax=Axes3D(fig)
ax.plot(x,y,z,'g',linewidth=2)
ax.set_zlim3d(-1.0,1.0)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_title('A spiral as a parametric curve',
         weight='bold',size=16)
ax.elev, ax.azim = 60, -120
fig.savefig('pic/fig-6-7-1.png')
'''

## 6.7.2. Lissajous曲线, 使用 mlab
# =============================================================================
# import numpy as np
# 
# theta=np.linspace(0,2*np.pi,101)
# a=0.3; m=11; n=9
# x=(1+a*np.cos(n*theta))*np.cos(m*theta)
# y=(1+a*np.cos(n*theta))*np.sin(m*theta)
# z=a*np.sin(n*theta)
# 
# from mayavi import mlab
# mlab.plot3d(x,y,z,np.sin(n*theta),
#             tube_radius=0.025,colormap='spectral')
# mlab.axes(line_width=2,nb_labels=5)
# mlab.title('A spiral wrapped around a torus',size=1.2)
# 
# =============================================================================



## 6.8.1. 使用mplot3d 可视化简单曲面
'''
import numpy as np
xx,yy=np.mgrid[-2:2:81j,-3:3:91j]
zz=np.exp(-2*xx**2-yy**2)*np.cos(2*xx)*np.cos(3*yy)

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

plt.ion()
fig=plt.figure()
ax=Axes3D(fig)
ax.plot_surface(xx,yy,zz,rstride=4,cstride=3,
                color='c',alpha=0.9)
ax.contour(xx,yy,zz,zdir='x',offset=-3.0,color='black')
ax.contour(xx,yy,zz,zdir='y',offset=4.0,color='blue')
ax.contour(xx,yy,zz,zdir='z',offset=-2.0)
ax.set_xlim3d(-3.0,2.0)
ax.set_ylim3d(-3.0,4.0)
ax.set_zlim3d(-2.0,1.0)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
#fig.savefig('pic/fig-6-8-1.png')
'''

## 6.8.2. 使用 mlab 可视化简单曲面

# =============================================================================
# import numpy as np
# 
# xx,yy=np.mgrid[-2:2:81j,-3:3:91j]
# zz=np.exp(-2*xx**2-yy**2)*np.cos(2*xx)*np.cos(3*yy)
# 
# from mayavi import mlab
# fig=mlab.figure()
# s=mlab.surf(xx,yy,zz,representation='surface')
# ax=mlab.axes(line_width=2,nb_labels=5)
# #mlab.title('Simple surface plot',size=0.4)
# 
# =============================================================================


## 6.9.1. Enneper 曲面，使用 mplot3d
#'''
import numpy as np
[u,v]=np.mgrid[-2:2:51j,-2:2:61j]
x,y,z=u*(1-u**2/3+v**2),v*(1-v**2/2+u**2),u**2-v**2

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

fig=plt.figure()
ax=Axes3D(fig)
ax.plot_surface(x.T,y.T,z.T,rstride=2,cstride=2,
    color='blue',alpha=0.5,linewidth=0.5)

ax.elev,ax.azim=50,-80
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_title('Enneper surface plot',
             weight='bold',size=14)
fig.savefig('pic/fig-6-9-1.png') 
#'''


## 6.9.2. Enneper 曲面，使用 mlab
# =============================================================================
# import numpy as np
# [u,v]=np.mgrid[-2:2:51j,-2:2:61j]
# x,y,z=u*(1-u**2/3+v**2),v*(1-v**2/2+u**2),u**2-v**2
# 
# from mayavi import mlab
# fig=mlab.figure()
# s=mlab.mesh(x,y,z,representation='surface',
#             linewidth=0.5,opacity=0.5)
# ax=mlab.axes(line_width=2,nb_labels=5)
# 
# =============================================================================


# =============================================================================
##6.10.  Julia Set, 使用 mlab, 参考 mayavi 官网
#'''
#import numpy as np
##set up initial grid
#x,y=np.ogrid[-1.5:0.5:1000j,-1.0:1.0:1000j]
#z=x+1j*y
#julia=np.zeros(z.shape)
#c=-0.7-0.4j
#
##build the Julia set
#for it in range(1,101):
#    z=z**2+c
#    escape=z*z.conj()>4
#    julia+=(1/float(it))*escape
#    
#from mayavi import mlab
#mlab.figure(size=(800,600))
#mlab.surf(julia,colormap='gist_ncar',
#          warp_scale='auto',vmax=1.5)
#mlab.view(15,30,500,[-0.5,-0.5,2.0])
#mlab.show()
# =============================================================================


## 6.10. Julia Set, 使用平面图像
'''
import numpy as np

x,y=np.mgrid[-1.5:1.0:1000j,-1.0:1.0:1000j]
z=x+1j*y
julia=np.zeros(z.shape)
c=-0.7-0.4j

for it in range(1,101):
    z=z*z+c
    escape=z*z.conjugate()>4
    julia+=(1/float(it))*escape

import matplotlib.pyplot as plt
plt.imshow(julia)
'''


  